Dans le cadre du module de statistiques appliquées, nous avions abordé les bases des statistiques et pris en main des outils permettant d’effectuer des opérations statistiques avancées grâce au langage R. Ainsi, notre choix s’est porté sur des données statistiques relatives au prestige (score attribuée) de différentes professions canadiennes.
Ce projet de statistiques appliquées consiste en l’utilisation d’outils R adaptés afin d’effectuer une modélisation statistique sur notre jeu de données et produire un rapport.
R est un langage de programmation dont le but est de pouvoir traiter et organiser des jeux de données afin de pouvoir y appliquer des tests statistiques plus ou moins complexes et se représenter ces données graphiquement à l’aide d’une grande variété de graphiques disponibles. RStudio est une application proposant un environnement de développement et des outils adaptés au langage et à l’environnement de programmation R.
La fonction principale de RStudio consiste à faciliter le développement d’applications en langage R. Pour ce faire, le programme dispose de nombreux outils qui vous permettent notamment de créer des scripts, compiler du code, créer des graphes, ainsi que de travailler avec divers jeux de données.
L’extension R markdown permet de générer des documents de manière dynamique en mélangeant texte mis en forme et résultats produits par du code R. Les documents générés peuvent être au format HTML, PDF, Word, et bien d’autres. C’est donc un outil très pratique pour l’exportation, la communication et la diffusion de résultats d’analyse.
Notre jeu de données comporte 102 individus décrits par 6 variables :
head(data) %>% paged_table()
viz_edu <- data %>%
ggplot(aes(x = prestige, y = education, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(6, 16, 2)
)
viz_edu <- ggplotly(viz_edu)
viz_inc <- data %>%
ggplot(aes(x = prestige, y = income, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(0, 20000, 10000)
)
viz_inc <- ggplotly(viz_inc)
viz_women <- data %>%
ggplot(aes(x = prestige, y = women, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(0, 100, 25)
)
viz_women <- ggplotly(viz_women)
viz_census <- data %>%
ggplot(aes(x = prestige, y = census, col = type)) +
geom_point() + theme_bw() +
scale_x_continuous(
breaks = seq(25, 75, 25)
) +
scale_y_continuous(
breaks = seq(2500, 7500, 2500)
)
viz_census <- ggplotly(viz_census)
subplot(viz_edu, viz_inc, viz_women, viz_census, nrows = 2, margin = 0.07)
En statistiques et en apprentissage automatique, un modèle de régression linéaire est un modèle de régression qui cherche à établir une relation linéaire entre une variable, dite expliquée, et une ou plusieurs variables, dites explicatives.
Dans un premier temps, nous pouvons étudier rapidement d’éventuelles corrélations entre deux variables à l’aide d’un nuage de points. Ici, nous utilisons la fonction ggpairs afin de pouvoir en observer plusieurs en même temps.
generate_scatterM <- function(data, mapping){
data %>% ggplot(mapping = mapping) + geom_point(size = 0.8) +
geom_smooth(method = loess, color = "red", size = 0.85, se = FALSE) +
geom_smooth(method = lm, color = "blue", size = 0.85, se = FALSE)
}
scatterM <- data %>% ggpairs(columns = 1:5, lower = list(continuous = generate_scatterM))
scatterM <- ggplotly(scatterM)
scatterM
Avec les observations précédentes, nous pouvons déjà écarter certaines linéarités et choisir lesquelles sont intéressantes à étudier. Ici, nous reprenons donc la relation entre la variable prestige (il s’agit d’un score de prestige relatif à la profession) et la variable education (qui reflète le niveau d’étude). Cette fois nous utilisons ggplot pour une meilleure représentation.
education <- data %>%
ggplot(aes(x = education, y = prestige)) +
geom_point() +
geom_smooth(method = loess, se = T) +
geom_smooth(method = lm, color = "red", se = F) +
theme_bw() +
scale_x_continuous(
breaks = seq(6, 16, 2)
) +
scale_y_continuous(
breaks = seq(20, 80, 20)
)
scatter_edu <- ggplotly(education)
scatter_edu
La droite bleue est définie par la méthode des moindres carrés (MSE), il s’agit d’une droite de régression linéaire entre les variables education et prestige.
La courbe rouge indique la tendance globale entre ces deux variables, il s’agit d’une courbe de régression de type lowess. Les deux courbes extérieures représentent l’intervalle de confiance de cette courbe de régression.
On constate que que la droite de régression est presque toujours comprise dans l’intervalle de confiance, l’hypothèse de linéarite entre les variables education et prestige est donc acceptable.
Nous testons aussi la relation entre les variables income et prestige.
income <- data %>%
ggplot(aes(x = income, y = prestige)) +
geom_point() +
geom_smooth(method = loess, se = T) +
geom_smooth(method = lm, color = "red", se = F) +
theme_bw() +
scale_x_continuous(
breaks = seq(0, 25000, 5000)
) +
scale_y_continuous(
breaks = seq(20, 100, 20)
)
scatter_inc <- ggplotly(income)
scatter_inc
Ici, en regardant la forme du lien entre la variable entre les deux variables, on s’aperçoit que la droite de régression suis bien moins l’intervalle de confiance de la courbe lowess. l’hypothèse de linéarité est alors plus critiquable.
Pour déterminer la droite de régression, on ajuste un modèle linéaire simple aux données, à l’aide de la fonction “lm”.
lin_model_edu <- lm(prestige ~ education, data = data)
lin_model_edu
##
## Call:
## lm(formula = prestige ~ education, data = data)
##
## Coefficients:
## (Intercept) education
## -10.732 5.361
« Intercept » correspond ici à l’ordonnée à l’origine, le « b » de notre droite, et le « x » est la pente de la droite, ce qui correspond au « b » dans notre notation.
L’équation, de notre droite est donc \(y = -10.732 + 5.361x\).
bacf <- acf(residuals(lin_model_edu), plot = F)
bacfdf <- with(bacf, data.frame(lag, acf))
conf.level <- 0.95
ciline <- qnorm((1 - conf.level)/2)/sqrt(length(residuals(lin_model_edu)))
lag_plot <- bacfdf %>%
ggplot(aes(x = lag, y = acf)) +
geom_hline(aes(yintercept = 0)) +
geom_segment(aes(xend = lag, yend = 0)) +
geom_hline(aes(yintercept = ciline), linetype = 3, color = 'darkblue') +
geom_hline(aes(yintercept = -ciline), linetype = 3, color = 'darkblue') +
theme_bw()
lag_plot <- ggplotly(lag_plot, tooltip = c("lag", "acf"))
lag_plot
durbinWatsonTest(lin_model_edu)
## lag Autocorrelation D-W Statistic p-value
## 1 0.2752512 1.43917 0.002
## Alternative hypothesis: rho != 0
qq_edu <- lin_model_edu %>%
ggplot(aes(sample = rstandard(.))) +
stat_qq_line(color = "red", linetype = "dashed") +
stat_qq(size = 1.2) +
theme_bw() +
scale_x_continuous(
breaks = seq(-2, 2, 1)
) +
scale_y_continuous(
breaks = seq(-3, 3, 1)
) +
labs(title = "Normal Q-Q",
x = "Theoretical Quantiles\nlm(prestige ~ education)",
y = "Standardized residuals")
qq_edu <- ggplotly(qq_edu)
qq_edu
shapiro.test(residuals(lin_model_edu))
##
## Shapiro-Wilk normality test
##
## data: residuals(lin_model_edu)
## W = 0.98065, p-value = 0.1406
lin_model_inc <- lm(prestige~income, data = data)
lin_model_inc
##
## Call:
## lm(formula = prestige ~ income, data = data)
##
## Coefficients:
## (Intercept) income
## 27.141176 0.002897
qq_inc <- lin_model_inc %>%
ggplot(aes(sample = rstandard(.)))+
stat_qq_line(color = "red", linetype = "dashed") +
stat_qq(size = 1.2) +
theme_bw() +
scale_x_continuous(
breaks = seq(-2, 2, 1)
) +
scale_y_continuous(
breaks = seq(-3, 3, 1)
) +
labs(title = "Normal Q-Q",
x = "Theoretical Quantiles\nlm(prestige ~ income)",
y = "Standardized residuals")
qq_inc <- ggplotly(qq_inc)
qq_inc
shapiro.test(residuals(lin_model_inc))
##
## Shapiro-Wilk normality test
##
## data: residuals(lin_model_inc)
## W = 0.97169, p-value = 0.02729
sc_loc_edu <- lin_model_edu %>%
ggplot(aes(fitted(.), sqrt(abs(rstandard(.))))) +
geom_point(size = 1.2) +
geom_smooth(method = loess, se = FALSE) +
theme_bw() +
scale_x_continuous(
breaks = seq(30, 70, 10)
) +
scale_y_continuous(
limits = c(0, 1.5),
breaks = seq(0, 1.5, 0.5)
) +
labs(title = "Scale-location",
x = "Fitted values\nlm(prestige ~ education)",
y = "sqrt |Standardized residuals|")
sc_loc_edu <- ggplotly(sc_loc_edu)
sc_loc_edu
ncvTest(lin_model_edu)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 0.6327545, Df = 1, p = 0.42635
res_v_fit_edu <- lin_model_edu %>%
ggplot(aes(fitted(.), residuals(.))) +
geom_point() +
stat_smooth(method="loess", se = FALSE) +
geom_hline(yintercept = 0, col = "red", linetype = "dashed") +
theme_bw() +
scale_x_continuous(
breaks = seq(30, 70, 10)
) +
scale_y_continuous(
limits = c(-30, 20),
breaks = seq(-30, 20, 10)
) +
labs(title = "Residuals vs. Fitted values",
x = "Fitted values\nlm(prestige ~ education)",
y = "Residuals")
res_v_fit_edu <- ggplotly(res_v_fit_edu)
res_v_fit_edu
summary(lin_model_edu)
##
## Call:
## lm(formula = prestige ~ education, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.0397 -6.5228 0.6611 6.7430 18.1636
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -10.732 3.677 -2.919 0.00434 **
## education 5.361 0.332 16.148 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.103 on 100 degrees of freedom
## Multiple R-squared: 0.7228, Adjusted R-squared: 0.72
## F-statistic: 260.8 on 1 and 100 DF, p-value: < 2.2e-16
confint(lin_model_edu) %>% kable()
| 2.5 % | 97.5 % | |
|---|---|---|
| (Intercept) | -18.027220 | -3.436744 |
| education | 4.702223 | 6.019533 |
predictions <- tibble(education=c(9.21, 11.07, 14.62))
predict(lin_model_edu, newdata = predictions, interval = "confidence") %>% kable()
| fit | lwr | upr |
|---|---|---|
| 38.64170 | 36.58966 | 40.69374 |
| 48.61293 | 46.81134 | 50.41452 |
| 67.64405 | 64.52387 | 70.76423 |
data$residuals <- residuals(lin_model_edu)
head(data) %>% kable()
| education | income | women | prestige | census | type | residuals |
|---|---|---|---|---|---|---|
| 13.11 | 12351 | 11.16 | 68.8 | 1113 | prof | 9.250875 |
| 12.26 | 25879 | 4.02 | 69.1 | 1130 | prof | 14.107621 |
| 12.77 | 9271 | 15.70 | 63.4 | 1171 | prof | 5.673573 |
| 11.42 | 8865 | 9.11 | 56.8 | 1175 | prof | 6.310758 |
| 14.62 | 8403 | 11.68 | 73.5 | 2111 | prof | 5.855950 |
| 15.64 | 11030 | 5.13 | 77.6 | 2113 | prof | 4.487854 |
data$fitted <- fitted(lin_model_edu)
head(data) %>% kable()
| education | income | women | prestige | census | type | residuals | fitted |
|---|---|---|---|---|---|---|---|
| 13.11 | 12351 | 11.16 | 68.8 | 1113 | prof | 9.250875 | 59.54913 |
| 12.26 | 25879 | 4.02 | 69.1 | 1130 | prof | 14.107621 | 54.99238 |
| 12.77 | 9271 | 15.70 | 63.4 | 1171 | prof | 5.673573 | 57.72643 |
| 11.42 | 8865 | 9.11 | 56.8 | 1175 | prof | 6.310758 | 50.48924 |
| 14.62 | 8403 | 11.68 | 73.5 | 2111 | prof | 5.855950 | 67.64405 |
| 15.64 | 11030 | 5.13 | 77.6 | 2113 | prof | 4.487854 | 73.11215 |
vcov(lin_model_edu) %>% kable()
| (Intercept) | education | |
|---|---|---|
| (Intercept) | 13.520978 | -1.1835055 |
| education | -1.183506 | 0.1102162 |
pred_interval <- predict(lin_model_edu, interval = "prediction")
data <- cbind(data, pred_interval)
linear_regression <- data %>%
ggplot(aes(y = prestige, x = education)) +
geom_point(size = 1.2) +
geom_smooth(color = "red", method = "lm", fill = "red") +
geom_line(aes(y = lwr), color = "red", linetype = "dashed") +
geom_line(aes(y = upr), color = "red", linetype = "dashed") +
annotate("text", x = 9, y = 80,
label = paste0("prestige = ", coef(lin_model_edu)["(Intercept)"] %>% round(3),
" + ", coef(lin_model_edu)["education"] %>% round(3), " * education")) +
theme_bw() +
scale_x_continuous(
breaks = seq(6, 16, 2)
) +
scale_y_continuous(
breaks = seq(25, 75, 25)
) +
labs(title = "Linear Regression",
x = "Education", y = "Prestige")
linear_regression <- ggplotly(linear_regression)
linear_regression